Данные

Данные об университетах.

Признаки:

Выберем некоторые из этих признаков.

uni <- read.csv2("/Users/Vladimir/Dropbox/Statistics_Autmn_2016/PCA/I.csv")
uni <- uni[,c("X", "PPIND", "NEW10", "R_B_COST","ADD_FEE", "BOOK", "SAL_ALL", "NUM_ALL" ,"APP_ACC", "SF_RATIO", "IN_STATE", "OUT_STAT")]
uni <- na.omit(uni)
nums <- c("NEW10", "R_B_COST","ADD_FEE", "BOOK", "SAL_ALL", "NUM_ALL" ,"APP_ACC", "SF_RATIO", "IN_STATE", "OUT_STAT")
uni[,nums] <- scale(log(uni[,nums]))

Матриксплот:

pairs.panels(uni[,nums],
               pch=21,
               lm=TRUE, 
               lwd = 1,
               ellipses = FALSE)

По скаттерплотам можно предположить, что существует два кластера. Это хорошо видно, например, на графике NUM_ALL vs. IN_STATE.

Model-based Clustering

Модель

Модель – смесь нормальных распределений: \[ p(x) = \sum_{k = 1}^G w_k p_k(x), ~p_k(x) = \phi(x; \mu_k, \Sigma_k). \]

Хотим найти параметры максимизируя функцию правдоподобия: \[ \mathcal{L}(\theta; \mathbf{x}) = \prod_{i = 1}^n \sum_{k = 1}^G w_k \phi(x_i; \mu_k, \Sigma_k). \]

Подробно разберем случай смеси двух распределений.

Смесь: \[ p(x) = (1 - \pi)\phi(x; \mu_1, \Sigma_1) + \pi\phi(x; \mu_2, \Sigma_2). \]

Логарифм функции правдоподобия:

\[ \mathcal{l}(\theta; \mathbf{x}) = \sum_{i = 1}^n log((1 - \pi)\phi(x_i; \mu_1, \Sigma_1) + \pi\phi(x_i; \mu_2, \Sigma_2)). \]

Численно максимизировать \(\mathcal{l}(\theta; x)\) сложно, поэтому вводятся ненаблюдаемые переменные \(\Delta_i\), принимающие значения \(0\) или \(1\). Если \(\Delta_i = 1\), то элемент выборки \(x_i\) пришел к нам из первого распределения, иначе – из второго. Логарифм функции правдопадобия перепишется так:

\[ \mathcal{l}(\theta; \mathbf{x}, \mathbf{\Delta}) = \sum_{i = 1}^n (1 - \Delta_i) \log \phi(x_i; \mu_1, \Sigma_1) + \Delta_i \log \phi(x_i; \mu_2, \Sigma_2) + \sum_{i = 1}^n (1 - \Delta_i) \log(1 - \pi) + \Delta_i \log \pi. \]

Так как значения \(\Delta_i\) неизвестны, заменим их на условные мат. ожидания:

\[ \gamma_i(\theta) = \mathbb{E}(\Delta_i|\theta, \mathbf{x}). \]

Далее воспользуемся ЕМ-алгоритмом, чтобы найти оценки параметров.

EM-алгоритм

Пусть у нас есть начальные оценки параметров: \(\hat{\mu}_1\), \(\widehat{\Sigma}_1\), \(\hat{\mu}_2\), \(\widehat{\Sigma}_2\), \(\hat{\pi}\).

Expectation Step

\[ \hat{\gamma}_i = \frac{\hat{\pi}\phi(x_i; \mu_2, \Sigma_2)}{(1 - \pi)\phi(x_i; \mu_1, \Sigma_1) + \pi\phi(x_i; \mu_2, \Sigma_2)}. \]

Maximization Step

\[ \hat{\mu}_1 = \frac{\sum_{i = 1}^n (1 - \hat{\gamma}_i)x_i}{\sum_{i = 1}^n (1 - \hat{\gamma}_i)}, \] \[ \hat{\mu}_2 = \frac{\sum_{i = 1}^n \hat{\gamma}_i x_i}{\sum_{i = 1}^n \hat{\gamma}_i}. \] Если не делать никаких предположений о виде ковариационных матриц, то их оценки: \[ \widehat{\Sigma}_1 = \frac{\sum_{i = 1}^n (1 - \hat{\gamma}_i) (x_i - \hat{\mu}_1)(x_i - \hat{\mu}_1)^{\mathrm{T}}}{\sum_{i = 1}^n (1 - \hat{\gamma}_i)}, \] \[ \widehat{\Sigma}_2 = \frac{\sum_{i = 1}^n \hat{\gamma}_i (x_i - \hat{\mu}_1)(x_i - \hat{\mu}_1)^{\mathrm{T}}}{\sum_{i = 1}^n \hat{\gamma}_i}, \] \[ \hat{\pi} = \frac{1}{n}\sum_{i = 1}^n \hat{\gamma}_i. \]

Продолжаем пока не сойдемся.

Ограничения на ковариационные матрицы

Ковариационные матрицы определяют вид эллипсов.

Найдем спектральное разложение ковариационной матрицы: \[ \Sigma_k = \lambda_k D_k A_k D_k^\mathrm{T}, \] где \(D_k\) – матрица, составленная из собственных векторов \(\Sigma_k\), \(A_k\) – диагональная матрица, компоненты которой пропорциональны собственным числам, \(\lambda_k\) – коэффициент пропорциональности.

Такое разложение можно интерпретировать геометрически: \(D_k\) отвечает за ориентацию эллипса, \(A_k\) – за его форму, \(\lambda_k\) – за его объем. Фиксируя различные признаки эллипсов, получаем различные модели. В случае размерности большей 2 модель обозначают тремя буквами, например, EVI означает, что у эллипсов одинаковый объем, возможно разные формы, а корреляции равны 0.

BIC

Байесовский информационный критерий.

Идея: давайте штрафовать за большое число параметров.

\[ BIC = -2 \hat{\mathcal{l}} + k \log(n), \] где \(n\) – число точек, \(k\) – число параметров

Кластеризация

raw.data <- uni[,nums]

Посмотрим на байесовский информационный критерий.

BIC <- mclustBIC(raw.data)
plot(BIC)

summary(BIC)
## Best BIC values:
##              EEE,2        EEE,5       EEE,3
## BIC      -2899.265 -2908.343604 -2913.00935
## BIC diff     0.000    -9.078216   -13.74397

BIC рекомендует предположить, что у нас два кластера, а ковариационные матрицы одинаковы. Но по матриксплотам выше видно, что точнее будет взять VEV, так как объемы разные, формы одинаковые, корреляции тоже разные (например, NEW_10 vs. R_B_COST).

m <- Mclust(raw.data, G = 2, modelNames = "VEV")

Результат классификации.

plot(m, "classification")

Картинка в главных компонентах

clusplot(raw.data, m$classification, main='2D representation of the Cluster solution',
         color=TRUE, shade=TRUE,
         labels=2, lines=0)

Silhouette plot

Идея – для каждой точки определяем насколько она похожа на остальные точки из своего кастера по сравнению с точками из других кластеров. Например, в смысле евклидовой метрики.

Пусть \(a(i)\) – среднее расстояние от \(i\)-й точки до остальных в кластере (cohesion), \(b(i)\) – самое маленькое из средних расстояний до точек из других кластеров (separation).

Силуэт определяется так:

\[ s(i) = \frac{b(i) - a(i)}{max\{a(i), b(i)\}}. \]

Чем ближе \(s(i)\) к \(1\), тем более уверенно можно сказать, что точка в нужном кластере. Среднее значение \(s\) – мера того, как хорошо кластеризовали.

dissE <- daisy(raw.data) 
dE2   <- dissE^2
sk2   <- silhouette(m$classification, dE2)
plot(sk2)

Кроме того, силуэт можно использовать как один из способов предположить количество кластеров. Если выбрали слишком большое или слишком маленькое k, некоторые селуэты окажутся заметно уже чем остальные. Например, если возьмем 6 кластеров для наших данных:

m6 <- Mclust(raw.data, G = 6, modelNames = "EEE")
dissE <- daisy(raw.data) #считаем среднее расстояние от каждой точки до остальных (внутри и вне кластера)
dE2   <- dissE^2
sk2   <- silhouette(m6$classification, dE2)
plot(sk2)

Partitioning Clustering

k-means

cl <- kmeans(raw.data, 2)
plot(raw.data, col = c("red", "blue")[cl$cluster])

Картинка в главных компонентах:

cl <- kmeans(scale(raw.data), 2)
clusplot(scale(raw.data), cl$cluster, main='2D representation of the Cluster solution',
         color=TRUE, shade=TRUE,
         labels=2, lines=0)

Silhouette plot

dissE <- daisy(raw.data) 
dE2   <- dissE^2
sk2   <- silhouette(cl$cluster, dE2)
plot(sk2)

Если верить силуэту, данные разделились чуть хуже чем при model-based clustering (среднее s(i) 0.5 против 0.51 в случае model-based).

Hierarchical Clustering

Complete Linkage

На каждом шаге считается расстояние до “дальнего соседа”.

\[ R(W, S) = \max_{s \in S, w \in W} \rho(w,s) \]

#euclid dist
d <- dist(raw.data)

hc.cl <- hclust(d, method = "complete")
plot(hc.cl)

plot(raw.data, col = c("red", "blue")[cutree(hc.cl, 2)])

Ward

Здесь расстояние – то, насколько увеличится сумма квадратов расстояний до центра при слиянии двух кластеров, цена слияния.

\[ \Delta(A, B) = \sum_{i \in A \cup B}\|x_i - \mu_{A \cup B} \|^2 - \sum_{i \in A}\|x_i - \mu_{A} \|^2 - \sum_{i \in B}\|x_i - \mu_{B} \|^2 \]

#euclid dist
d <- dist(raw.data)

hc.ward <- hclust(d, method = "ward.D2")
plot(hc.ward)

Так как метод Варда на каждом шаге выбирает минимальную цену слияния, большие скачки могут означать, что мы объединяем точки не внутри кластеров, а сами кластеры. Тогда по картинке выше можно решить, что у нас два кластера.

plot(raw.data, col = c("red", "blue")[cutree(hc.ward, 2)])

Проверка

На самом деле университеты делятся на государственные и частные. Выглядит это так:

pairs.panels(uni[,nums],
               pch=21 +as.numeric(uni$PPIND),
               bg = c("red", "blue")[uni$PPIND],
               lm=TRUE, 
               lwd = 1,
               ellipses = FALSE)

Проверим, насколько ошиблись методы выше.

Model-Based

table(uni$PPIND, m$classification)
##    
##      1  2
##   1 75  0
##   2  3 42

Всего три ошибки.

k-means

table(uni$PPIND, cl$cluster)
##    
##      1  2
##   1  2 73
##   2 42  3

Ошибок чуть больше чем в model-based. Если вернуться к матриксплоту выше, видно, что model-based справился, а k-means – не совсем.

Hierarchical Clustering

Complete linkage

table(uni$PPIND, cutree(hc.cl, 2))
##    
##      1  2
##   1 67  8
##   2  2 43

Ward

table(uni$PPIND, cutree(hc.cl, 2))
##    
##      1  2
##   1 67  8
##   2  2 43

Количество ошибок больше чем в двух методах выше. По матриксплотам выше видно, что в случае иерархической кластеризации получилось что-то похожее на k-means.